1.1 用TTree 结构存储事件数据

目的:学习如何将事件以TTree的格式存储到ROOT文件

内容:

  1. 束流入射到靶上产生核反应,靶上生成的中子和$\gamma$入射到长条型塑料闪烁体探测器的不同位置,并在不同深度与探测介质产生相互作用。
  2. 入射粒子形成的光信号向探测器两端传输,形成时间信号$t_u$和$t_d$, 能量信号 $Q_u, Q_d$。
  3. 将入射粒子的原始信息,以及探测器信息逐事件存入root文件的TTree结构中.
    • 入射粒子原始信息:粒子种类、能量,粒子在探测器上的入射位置、入射深度等
    • 探测器信息: 两端时间信息$t_u$、$t_d$,能量沉积信息$q_u$,$q_d$
  4. 打开root文件,进行基本的分析。

中子飞行时间谱测量原理

探测器参数

探测器的几何设置见下图:

  • 探测器长度:2$L$, 忽略探测器宽度。
  • 探测器厚度:$\Delta D$
  • 光在闪烁体的传播速度: $v_{sc}$
  • 光在闪烁体的衰减常数:$\lambda$
  • 粒子从起始位置到探测器中心的距离:$D$
  • 粒子入射位置:$x$
  • 粒子相互作用的位置(深度):$(-\frac{\Delta D}{2} , +\frac{\Delta D}{2})$

公式

飞行距离:

$$ d=\sqrt{D^2+x^2} $$

飞行时间, 以靶为起点:

$$ \begin{equation} 中子:TOF_n(ns)=\frac{72.29824 \cdot d(m)}{\sqrt{E_n(MeV)}}\\ \gamma:TOF_{\gamma}(ns)=3.333\cdot d(m) \end{equation} $$

探测器两端时间信号(不考虑由电缆长度或PMT度越时间引起的offset):

$$ \begin{equation} t_u=TOF+(L-x)/v_{sc}\\ t_d=TOF+(L+x)/v_{sc}\\ x=(t_d-t_u)*v_{sc}/2\\ TOF=(t_u+t_d)/2 -L/v_{sc}\\ \end{equation} $$

探测器两端能量,设x 处能量沉积 $2Q_0$:

$$ Q_u=Q_0 exp[-(L-x)/\lambda]\\ Q_d=Q_0 exp[-(L+x)/\lambda]\\ x=\frac{\lambda}{2}ln\frac{q_u}{q_d}\\ Q_0=e^{2L/\lambda}*Q_u*Q_d\\ $$

探测器分辨率:

  • 时间分辨率(FWHM):$T_{res.}$
  • 相对能量分辨率$\Delta E/E$: $E_{res.}$

随机数:

  • TRandom3 class
    • Uniform(Double_t x1, Double_t x2);
    • Gaus(Double_t mean, Double_t sigma)

代码

// 运行方法:
// 将下列代码保存到tree.cc
// 在ROOT环境中运行:
// root -l //进入ROOT环境
// .x tree.cc //运行脚本
// .q //退出root环境
// root -l tree.root //打开ROOT文件
void tree(){ 
//常量声明
  const Double_t D=500.;//cm, distance between target and the scin.(Center)
  const Double_t L=100.;//cm, half length of the scin.
  const Double_t dD=5.;//cm, thickness of the scin.
  const Double_t TRes=1.;//ns, time resolution(FWHM) of the scintillator.
  const Double_t Lambda=380.;//cm, attenuation lenght of the scin.
  const Double_t QRes=0.1;//relative energy resolution(FWHM) of the scin. 
  const Double_t Vsc=7.5;//ns/cm, speed of light in the scin.
  const Double_t En0=100;//MeV, average neutron energy
  const Double_t EnRes=50.;//MeV, energy spread of neutron(FWHM)
  const Double_t Eg0=1;//MeV, gamma energy  
  const Double_t Rg=0.3;//ratio of gamma,ratio of neutron 1-Rg 

  //1. 声明在tree结构中定义需要的变量分支
  Double_t x;//入射位置
  Double_t e;//能量
  int pid;    //粒子种类,n:pid=1,g:pid=0
  Double_t tof, ctof;//TOF:粒子实际飞行时间,cTOF:计算得到的TOF
  Double_t tu, td;
  Double_t qu, qd;

  Double_t tu_off=5.5;//time offset,//PMT的度越时间+电缆上的传输时间
  Double_t td_off=20.4;//time offset 

 //2. 定义新ROOT文件,声明新的Tree 
  TFile *opf=new TFile("tree.root","recreate");//新文件tree.root,指针 *opf
  TTree *opt=new TTree("tree","tree structure");//新tree,指针 *opt

  //3. 将变量地址添加到tree结构中
    //第一个参数为变量名称,第二个为上面定义的变量地址,第三个为变量的类型说明,D表示Double_t,I表示 Int_t。
  opt->Branch("x", &x, "x/D");
  opt->Branch("e", &e, "e/D");
  opt->Branch("tof", &tof, "tof/D");
  opt->Branch("ctof",&ctof,"ctof/D");
  opt->Branch("pid", &pid, "pid/I");
  opt->Branch("tu", &tu, "tu/D");
  opt->Branch("td", &td, "td/D");
  opt->Branch("qu", &qu, "qu/D"); 
  opt->Branch("qd", &qd, "qd/D");  

// histogram,ROOT文件中除了TTree结构外,还可存储histogram,graph等
  TH1D *hctof=new TH1D("hctof","neutron time of flight",1000,0,100);
  TRandom3 *gr=new TRandom3(0);//声明随机数

  //4. 循环,计算变量的值,逐事件往tree结构添加变量值。
  for(int i=0;i<100000;i++){
    x=gr->Uniform(-L, L);//粒子入射位置,在-L,L范围内均匀抽样.
    Double_t Dr=D+gr->Uniform(-0.5,0.5)*dD;//粒子在探测器厚度范围内均匀产生光信号
    Double_t d=TMath::Sqrt(Dr*Dr+x*x);//粒子实际飞行距离,flight path
    if(gr->Uniform() < Rg) { //判断为gamma入射
       pid=0;
       e=Eg0;
       tof=3.333*(d*0.01);
    }
    else {  //neutron
        pid=1;
        e=gr->Gaus(En0, EnRes/2.35); // energy of neutrons
        tof=72.29824/TMath::Sqrt(e)*(d*0.01);//ns
    }
    tu=tof+(L-x)/Vsc+gr->Gaus(0,TRes/2.35)+tu_off;
    td=tof+(L+x)/Vsc+gr->Gaus(0,TRes/2.35)+td_off;
    ctof=(tu+td)/2.;//simplified calculation.
    hctof->Fill(ctof);
//neutron:energy of recoil proton in plas. q0=0-En; gamma:q0=0-Egamma,compton plateau
    Double_t q0=e*gr->Uniform();
    qu=q0*TMath::Exp(-(L-x)/Lambda);
    qu=gr->Gaus(qu,qu*QRes/2.35);//QRes, relative energy resolution
    qd=q0*TMath::Exp(-(L+x)/Lambda);
    qd=gr->Gaus(qd,qd*QRes/2.35);
    opt->Fill();//5.将计算好的变量值填到Tree中

    if(i%1000==0) cout<<".";
  }
  cout<<endl;
  // 6.将数据写入root文件中
  hctof->Write();//写入预定义的histogram到文件
  opt->Write();//写入TTree到文件
  opf->Close();//关闭文件
}

练习:

  • 理解上述1.-6.步骤的逻辑关系。自行练习如何将数据已TTree格式存到ROOT文件中。
  • 理解如何在模拟中加入探测器分辨率信息。

查看ROOT文件

  • juypter运行tree->Draw()时,需要有网络链接才能正常显示。
In [1]:
//%jsroot on
In [2]:
TFile *ipf=new TFile("tree.root");//root -l tree.root
ipf->ls()//ROOT 环境下 > .ls
TFile**		tree.root	
 TFile*		tree.root	
  KEY: TH1D	hctof;1	neutron time of flight
  KEY: TTree	tree;1	tree structure
In [3]:
//观察tree的结构
tree->Print()
******************************************************************************
*Tree    :tree      : tree structure                                         *
*Entries :   100000 : Total =         6822821 bytes  File  Size =    5757076 *
*        :          : Tree compression factor =   1.18                       *
******************************************************************************
*Br    0 :x         : x/D                                                    *
*Entries :   100000 : Total  Size=     802585 bytes  File Size  =     614377 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.31     *
*............................................................................*
*Br    1 :e         : e/D                                                    *
*Entries :   100000 : Total  Size=     802585 bytes  File Size  =     565372 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.42     *
*............................................................................*
*Br    2 :tof       : tof/D                                                  *
*Entries :   100000 : Total  Size=     802645 bytes  File Size  =     742186 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.08     *
*............................................................................*
*Br    3 :ctof      : ctof/D                                                 *
*Entries :   100000 : Total  Size=     802675 bytes  File Size  =     741725 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.08     *
*............................................................................*
*Br    4 :pid       : pid/I                                                  *
*Entries :   100000 : Total  Size=     401467 bytes  File Size  =      42316 *
*Baskets :       13 : Basket Size=      32000 bytes  Compression=   9.47     *
*............................................................................*
*Br    5 :tu        : tu/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     755814 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.06     *
*............................................................................*
*Br    6 :td        : td/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     753729 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.06     *
*............................................................................*
*Br    7 :qu        : qu/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     769515 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.04     *
*............................................................................*
*Br    8 :qd        : qd/D                                                   *
*Entries :   100000 : Total  Size=     802615 bytes  File Size  =     769558 *
*Baskets :       26 : Basket Size=      32000 bytes  Compression=   1.04     *
*............................................................................*

观察指定事件数据

In [4]:
tree->Show(0);//显示指定的事件数据 event=1
======> EVENT:0
 x               = -50.0225
 e               = 95.2151
 tof             = 37.0316
 ctof            = 63.0591
 pid             = 1
 tu              = 62.6454
 td              = 63.4728
 qu              = 9.07882
 qd              = 10.9855
In [5]:
tree->Show(100);
======> EVENT:100
 x               = -91.9996
 e               = 112.608
 tof             = 34.5339
 ctof            = 61.1095
 pid             = 1
 tu              = 66.1766
 td              = 56.0423
 qu              = 49.7471
 qd              = 78.3448
In [6]:
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
c1->Clear();//* 在ROOT环境下可省略
tree->Draw("td-tu>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
In [7]:
tree->Draw("tof>>hh(1000,0,100)");//实际飞行时间
c1->SetLogy();
c1->Draw();
In [8]:
TH1D *hh=(TH1D*) ipf->Get("hctof"); //在代码中得到文件内数据指针的方法,在ROOT环境下课直接用 hctof->Draw()
hh->Draw();//显示文件内histogram
c1->Draw();

思考题1:观察上两个图有什么不同,分析其原因

In [9]:
c1->SetLogy(0);
gStyle->SetPalette(1);
tree->Draw("ctof:td-tu>>(2000,-20,50,1000,0,100)","","colz");//二维图显示
c1->Draw();

检查参数之间的关系是否与物理条件一致

  • gamma的到达探测器的时间(飞行时间)与探测器的x位置无关?
In [10]:
tree->Draw("ctof:td-tu>>(2000,-20,50,50,40,45)","pid==0","colz");//加入条件,gamma
c1->Draw();//为什么在两侧上弯?
In [11]:
tree->Draw("td-tu:x","","colz");//观察两个参量的关联
c1->Draw();
In [12]:
tree->Draw("log(qu/qd):x","","colz");
c1->Draw();
In [13]:
!jupyter nbconvert 1.1_create_tree.ipynb --to html
[NbConvertApp] Converting notebook 1.1_create_tree.ipynb to html

[NbConvertApp] Writing 442000 bytes to 1.1_create_tree.html

In [ ]: